/*---------------------------------------------------------------------------*\

    DAFoam  : Discrete Adjoint with OpenFOAM
    Version : v2

\*---------------------------------------------------------------------------*/

#include "checkGeometry.H"

namespace Foam
{
//- Default transformation behaviour for position
class transformPositionList
{
public:
    //- Transform patch-based field
    void operator()(
        const coupledPolyPatch& cpp,
        List<pointField>& pts) const
    {
        // Each element of pts is all the points in the face. Convert into
        // lists of size cpp to transform.

        List<pointField> newPts(pts.size());
        forAll(pts, facei)
        {
            newPts[facei].setSize(pts[facei].size());
        }

        label index = 0;
        while (true)
        {
            label n = 0;

            // Extract for every face the i'th position
            pointField ptsAtIndex(pts.size(), Zero);
            forAll(cpp, facei)
            {
                const pointField& facePts = pts[facei];
                if (facePts.size() > index)
                {
                    ptsAtIndex[facei] = facePts[index];
                    n++;
                }
            }

            if (n == 0)
            {
                break;
            }

            // Now ptsAtIndex will have for every face either zero or
            // the position of the i'th vertex. Transform.
            cpp.transformPosition(ptsAtIndex);

            // Extract back from ptsAtIndex into newPts
            forAll(cpp, facei)
            {
                pointField& facePts = newPts[facei];
                if (facePts.size() > index)
                {
                    facePts[index] = ptsAtIndex[facei];
                }
            }

            index++;
        }

        pts.transfer(newPts);
    }
};
}

bool Foam::checkCoupledPoints(
    const polyMesh& mesh,
    const bool report,
    labelHashSet* setPtr)
{
    const pointField& p = mesh.points();
    const faceList& fcs = mesh.faces();
    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    // Zero'th point on coupled faces
    //pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
    List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());

    // Exchange zero point
    forAll(patches, patchi)
    {
        if (patches[patchi].coupled())
        {
            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(
                patches[patchi]);

            forAll(cpp, i)
            {
                label bFacei = cpp.start() + i - mesh.nInternalFaces();
                const face& f = cpp[i];
                nbrPoints[bFacei].setSize(f.size());
                forAll(f, fp)
                {
                    const point& p0 = p[f[fp]];
                    nbrPoints[bFacei][fp] = p0;
                }
            }
        }
    }
    syncTools::syncBoundaryFaceList(
        mesh,
        nbrPoints,
        eqOp<pointField>(),
        transformPositionList());

    // Compare to local ones. Use same tolerance as for matching
    label nErrorFaces = 0;
    scalar avgMismatch = 0;
    label nCoupledPoints = 0;

    forAll(patches, patchi)
    {
        if (patches[patchi].coupled())
        {
            const coupledPolyPatch& cpp =
                refCast<const coupledPolyPatch>(patches[patchi]);

            if (cpp.owner())
            {
                scalarField smallDist(
                    cpp.calcFaceTol(
                        //cpp.matchTolerance(),
                        cpp,
                        cpp.points(),
                        cpp.faceCentres()));

                forAll(cpp, i)
                {
                    label bFacei = cpp.start() + i - mesh.nInternalFaces();
                    const face& f = cpp[i];

                    if (f.size() != nbrPoints[bFacei].size())
                    {
                        FatalErrorInFunction
                            << "Local face size : " << f.size()
                            << " does not equal neighbour face size : "
                            << nbrPoints[bFacei].size()
                            << abort(FatalError);
                    }

                    label fp = 0;
                    forAll(f, j)
                    {
                        const point& p0 = p[f[fp]];
                        scalar d = mag(p0 - nbrPoints[bFacei][j]);

                        if (d > smallDist[i])
                        {
                            if (setPtr)
                            {
                                setPtr->insert(cpp.start() + i);
                            }
                            nErrorFaces++;

                            break;
                        }

                        avgMismatch += d;
                        nCoupledPoints++;

                        fp = f.rcIndex(fp);
                    }
                }
            }
        }
    }

    reduce(nErrorFaces, sumOp<label>());
    reduce(avgMismatch, maxOp<scalar>());
    reduce(nCoupledPoints, sumOp<label>());

    if (nCoupledPoints > 0)
    {
        avgMismatch /= nCoupledPoints;
    }

    if (nErrorFaces > 0)
    {
        if (report)
        {
            Info << "  **Error in coupled point location: "
                 << nErrorFaces
                 << " faces have their 0th or consecutive vertex not opposite"
                 << " their coupled equivalent. Average mismatch "
                 << avgMismatch << "."
                 << endl;
        }

        return true;
    }
    else
    {
        if (report)
        {
            Info << "    Coupled point location match (average "
                 << avgMismatch << ") OK." << endl;
        }

        return false;
    }
}

Foam::label Foam::checkGeometry(
    const polyMesh& mesh,
    const autoPtr<surfaceWriter>& surfWriter,
    const autoPtr<writer<scalar>>& setWriter,
    const label maxIncorrectlyOrientedFaces)
{
    label noFailedChecks = 0;

    //Info << "\nChecking geometry..." << endl;

    // Get a small relative length from the bounding box
    const boundBox& globalBb = mesh.bounds();

    Info << "    Overall domain bounding box "
         << globalBb.min() << " " << globalBb.max() << endl;

    // Geometric directions
    const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one);
    labelList tmpList(3);
    forAll(tmpList, idxI) tmpList[idxI] = validDirs[idxI] / 2;
    Info << "    Mesh has " << mesh.nGeometricD()
         << " geometric (non-empty/wedge) directions " << tmpList << endl;

    // Solution directions
    const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one);
    forAll(tmpList, idxI) tmpList[idxI] = solDirs[idxI] / 2;
    Info << "    Mesh has " << mesh.nSolutionD()
         << " solution (non-empty) directions " << tmpList << endl;

    if (mesh.nGeometricD() < 3)
    {
        FatalErrorInFunction
            << "Mesh geometric directions is less than 3 and not supported!"
            << abort(FatalError);
    }

    if (mesh.checkClosedBoundary(true))
    {
        noFailedChecks++;
    }

    {
        cellSet cells(mesh, "nonClosedCells", mesh.nCells() / 100 + 1);
        cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells() / 100 + 1);
        if (
            mesh.checkClosedCells(
                true,
                &cells,
                &aspectCells,
                mesh.geometricD()))
        {
            noFailedChecks++;

            label nNonClosed = returnReduce(cells.size(), sumOp<label>());

            if (nNonClosed > 0)
            {
                Info << "  <<Writing " << nNonClosed
                     << " non closed cells to set " << cells.name() << endl;
                cells.instance() = mesh.pointsInstance();
                cells.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), cells);
                }
            }
        }

        label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());

        if (nHighAspect > 0)
        {
            Info << "  <<Writing " << nHighAspect
                 << " cells with high aspect ratio to set "
                 << aspectCells.name() << endl;
            aspectCells.instance() = mesh.pointsInstance();
            aspectCells.write();
            if (surfWriter.valid())
            {
                mergeAndWrite(surfWriter(), aspectCells);
            }
        }
    }

    {
        faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces() / 100 + 1);
        if (mesh.checkFaceAreas(true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info << "  <<Writing " << nFaces
                     << " zero area faces to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), faces);
                }
            }
        }
    }

    {
        cellSet cells(mesh, "zeroVolumeCells", mesh.nCells() / 100 + 1);
        if (mesh.checkCellVolumes(true, &cells))
        {
            noFailedChecks++;

            label nCells = returnReduce(cells.size(), sumOp<label>());

            if (nCells > 0)
            {
                Info << "  <<Writing " << nCells
                     << " zero volume cells to set " << cells.name() << endl;
                cells.instance() = mesh.pointsInstance();
                cells.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), cells);
                }
            }
        }
    }

    {
        faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces() / 100 + 1);
        if (mesh.checkFaceOrthogonality(true, &faces))
        {
            noFailedChecks++;
        }

        label nFaces = returnReduce(faces.size(), sumOp<label>());

        if (nFaces > 0)
        {
            Info << "  <<Writing " << nFaces
                 << " non-orthogonal faces to set " << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
            if (surfWriter.valid())
            {
                mergeAndWrite(surfWriter(), faces);
            }
        }
    }

    {
        faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces() / 100 + 1);
        if (mesh.checkFacePyramids(true, -SMALL, &faces))
        {
            if (maxIncorrectlyOrientedFaces > 0)
            {
                Info << "maxIncorrectlyOrientedFaces threshold is set to " << maxIncorrectlyOrientedFaces << endl;
            }

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            // if nFaces is less than maxIncorrectlyOrientedFaces
            // we do not consider it is a failed mesh. This allows
            // users to relax the mesh check for incorrectly oriented
            // mesh faces
            if (nFaces > maxIncorrectlyOrientedFaces)
            {
                noFailedChecks++;
            }

            if (nFaces > 0)
            {
                Info << "  <<Writing " << nFaces
                     << " faces with incorrect orientation to set "
                     << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), faces);
                }
            }
        }
    }

    {
        faceSet faces(mesh, "skewFaces", mesh.nFaces() / 100 + 1);
        if (mesh.checkFaceSkewness(true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info << "  <<Writing " << nFaces
                     << " skew faces to set " << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), faces);
                }
            }
        }
    }

    {
        faceSet faces(mesh, "coupledFaces", mesh.nFaces() / 100 + 1);
        if (checkCoupledPoints(mesh, true, &faces))
        {
            noFailedChecks++;

            label nFaces = returnReduce(faces.size(), sumOp<label>());

            if (nFaces > 0)
            {
                Info << "  <<Writing " << nFaces
                     << " faces with incorrectly matched 0th (or consecutive)"
                     << " vertex to set "
                     << faces.name() << endl;
                faces.instance() = mesh.pointsInstance();
                faces.write();
                if (surfWriter.valid())
                {
                    mergeAndWrite(surfWriter(), faces);
                }
            }
        }
    }

    return noFailedChecks;
}
